# Replication materials for Dür, Huber & Mateo Compensating the Losers: The (Limited) Elite-Public Gap in Trade Politics, *European Journal of Political Research*

# Load packages and set working directory ####

library(tidyverse)
library(labelled)
library(ordinal)
library(modelsummary)
library(ggstats)

here::i_am("EJPR_Replication_Materials.Rproj")

#The `marginaleffects` package changed its structure. This replication code works with version 0.23.0. The following chunk ensures the correct version is loaded from CRAN.
if (utils::packageVersion("marginaleffects") >= "0.23.0") {
  message("This code replicates with version 0.23.0 of marginaleffects but may not work with newer data")
} else { 
  
  remotes::install_version("marginaleffects", version = "0.23.0")
  
  message("This code is tested for version 0.23.0 of the marginaleffects package")}

# Dataset creation ####

#Here we load the data. The not-shareable age variable of legislators is stored in a separate datafile, which is not shared in this replication materials. Hence, you cannot reproduce the numerical results with the shared data for all analyses that require age. They are commented out here. The Analysis.html file provides a log with the analyses.


df0 <- readRDS("data/dat_legislators.rds") %>%
  mutate(ideology = ifelse(party_econ_combined == 0, 1,
                           ifelse(party_econ_combined == 6, 5,
                                  ifelse(is.na(party_econ_combined), NA, party_econ_combined))),
         university_education = ifelse(university_education == "Yes", "yes", university_education),
         gender = factor(ifelse(gender== "male", "Male", "Female")),
         q3_treat = as.factor(ifelse(q3_treat=="Businesses", "firms", "workers")))

if (dir.exists("not_shareable/")){
df0 <- dplyr::left_join(df0,
                        readRDS("not_shareable/dat_age.rds"),
                        by = "id")} else{
                          
                      message("To protect respondents privacy, the age data of legislators is excluded. Instead we add an empty age variable.")
                      df0$age <- NA}

df1 <- readRDS("data/dat_ES_PL_citizens.RDS") %>% 
  mutate(uni = ifelse(uni == "University", "yes", "no"),
         country = ifelse(country == "ES", "Spain", "Poland"),
         ideology = ifelse(lrecon <= 1, 1, 
                           ifelse(lrecon %in% c(2, 3), 2, 
                                  ifelse(lrecon %in% c(4, 5 , 6), 3,
                                         ifelse(lrecon %in% c(7, 8), 4, 5))))) %>% 
  select("country", "target", "subsidies", "tariffs", "tax", "infrastructure", 
         "gender", "age", "ideology", "uni", "employ")

df1 <- as.data.frame(unclass(df1))
df1 <- df1 %>%
  dplyr::mutate(across(c("subsidies", "tariffs", "tax", "infrastructure"), ~  . - 1))

df2 <- readRDS("data/data_GB_citizens.RDS") %>% 
  mutate(uni = ifelse(edu=="Bachelor or other university degree", "yes", "no"),
         country = "United Kingdom",
         ideology = ifelse(lrecon <= 1, 1, 
                           ifelse(lrecon %in% c(2, 3), 2, 
                                  ifelse(lrecon %in% c(4, 5 , 6), 3,
                                         ifelse(lrecon %in% c(7, 8), 4, 5))))) %>% 
  select("treatment", "subsidies", "tariffs", "tax", "infrastructure", 
         "country", "gender", "age", "ideology", "uni", "employ")

df2 <- as.data.frame(unclass(df2))
df2 <- df2 %>%
  dplyr::mutate(across(c("subsidies", "tariffs", "tax", "infrastructure"), ~  . - 1))


df <- data.frame(country=c(df0$country, df1$country, df2$country),
                 type=c(rep("legislator", nrow(df0)), rep("public", nrow(df1) + nrow(df2))),
                 subsidies = c(df0$q3_subsidies, df1$subsidies, df2$subsidies),
                 tariffs = c(df0$q3_tariffs, df1$tariffs, df2$tariffs),
                 tax = c(df0$q3_tax, df1$tax, df2$tax),
                 infrastructure = c(df0$q3_infrastructure, df1$infrastructure, df2$infrastructure),
                 age = c(df0$age, df1$age, df2$age),
                 gender = c(df0$gender, df1$gender, df2$gender),
                 university = c(df0$university_education, df1$uni, df2$uni),
                 treatment = c(df0$q3_treat, df1$target, df2$treatment), 
                 ideology =c(df0$ideology, df1$ideology, df2$ideology),
                 employ = factor(c(rep("Employed", nrow(df0)), as.character(df1$employ), as.character(df2$employ)),
                                 levels = c("Unemployed", "Employed")))


europe <- c("Austria", "Belgium", "Denmark", "Finland", "France", "Germany", "Hungary", "Iceland", 
            "Ireland", "Italy", "Luxembourg", "Moldova","Norway", "Portugal", "Romania",
            "Spain", "Sweden", "Switzerland", "United Kingdom", "Poland")

df3 <- df[df$country %in% europe, ]

# Descriptive figures ####

# .. Figure 1 ####


df3 %>% 
  filter(type=="legislator") %>%
  group_by(country) %>% 
  summarise(count= n()) %>% 
  ggplot(aes(x=reorder(country, -count), y=count)) +
  geom_bar(stat="identity", fill = "lightgray", colour = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 8)) +
  ylab("Count") +
  xlab("") +
  NULL

ggsave("figures/Figure1.pdf",
       width = 16, height = 9, units = "cm")

# .. Figure 2 ####


df3 %>% 
  select(type, subsidies, tariffs, tax, infrastructure) %>% 
  mutate(
    type = factor(type, levels = c("legislator", "public"), 
                  labels = c("Legislators", "Public")),
    subsidies = factor(subsidies, levels = 0:4, 
                       labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support")),
    tariffs = factor(tariffs, levels = 0:4, 
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support")),
    tax = factor(tax, levels = 0:4, 
                 labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support")),
    infrastructure = factor(infrastructure, levels = 0:4, 
                            labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))
  ) %>% 
  ggstats::gglikert(., c(5,2,3,4), 
                    facet_cols = vars(type), 
                    reverse_likert = FALSE, 
                    add_totals = TRUE,
                    totals_include_center = FALSE,
                    variable_labels=c(subsidies = "Subsidies",
                                      tariffs = "Tariffs",
                                      tax = "Tax cuts",
                                      infrastructure = "Infrastructure projects"),
                    labels_size = 2,
                    y_label_wrap = 10) +
  theme(legend.text = element_text(size = 8)) +
  scale_fill_brewer(palette = "RdYlBu")

ggsave("figures/figure2.pdf",
       width = 21.3333, height = 12, units = "cm")

# .. Figure 3 ####

df3 %>% 
  select(type, ideology, subsidies, tariffs, tax, infrastructure) %>% 
  mutate(
    type = factor(type, levels = c("legislator", "public"), 
                  labels = c("Legislators", "Public")),
    ideology = factor(ifelse(ideology < 3, "Left-wing",
                             ifelse(ideology > 3, "Right-wing", 
                                    ifelse(ideology == 3, "Centre", NA))),
                      levels = c("Left-wing", "Centre", "Right-wing")),
    subsidies = factor(subsidies, levels = 0:4, 
                       labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support")),
    tariffs = factor(tariffs, levels = 0:4, 
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support")),
    tax = factor(tax, levels = 0:4, 
                 labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support")),
    infrastructure = factor(infrastructure, levels = 0:4, 
                            labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))
  ) %>% 
  filter(!is.na(ideology)) %>% 
  ggstats::gglikert(., c(6,3, 4,5), 
                    facet_cols = vars(type), 
                    facet_rows = vars(ideology), 
                    reverse_likert = FALSE, 
                    add_totals = TRUE,
                    totals_include_center = FALSE,
                    variable_labels=c(subsidies = "Subsidies",
                                      tariffs = "Tariffs",
                                      tax = "Tax cuts",
                                      infrastructure = "Infrastructure projects"),
                    labels_size = 2,
                    y_label_wrap = 10) +
  theme(legend.text = element_text(size = 8)) +
  scale_fill_brewer(palette = "RdYlBu")

ggsave("figures/figure3.pdf",
       width = 21.3333, height = 16, units = "cm")

# .. Table A1 ####


#ideology
with(df3[df3$type=="public",], table(ideology))
with(df3[df3$type=="public",], prop.table(table(ideology)))
with(df3[df3$type=="legislator",], table(ideology))
with(df3[df3$type=="legislator",], prop.table(table(ideology)))
with(df3[df3$type=="legislator",], summary(ideology))

#gender
prop.table(table(df3$gender, df3$type), 2)

#education
prop.table(table(df3$university, df3$type), 2)

#treatment
prop.table(table(df3$treatment, df3$type), 2)

#all variables
df3$legislator <- ifelse(df3$type=="legislator", 1, 0)
df3$male <- ifelse(df3$gender=="Male", 1, 0)
df3$diverse <- ifelse(df3$gender=="Other", 1, 0)
df3$university_ed <- ifelse(df3$university=="yes", 1, 0)
df3$employed1 <- ifelse(df3$employ=="Employed", 1, 0)
df3$treatment_workers <- ifelse(df3$treatment=="workers", 1, 0)
df3$type_nice <- ifelse(df3$type=="legislator", "Legislator", "Public")

if (dir.exists("not_shareable/")) {
  datasummary((Heading("Infrastructure projects") * infrastructure)*(Heading("Actor type") * type_nice) +
              (Heading("Subsidies") * subsidies)*(Heading("Actor type") * type_nice) +
              (Heading("Tariff increases") * tariffs)*(Heading("Actor type") * type_nice) +
              (Heading("Tax cuts") * tax)*(Heading("Actor type") * type_nice) +
              #Heading("Legislator") * legislator)*(Heading("Actor type") * type) +
              (Heading("Political ideology") * ideology)*(Heading("Actor type") * type_nice) +
              (Heading("Age") * age)*(Heading("Actor type") * type_nice) + 
              (Heading("Male") * male)*(Heading("Actor type") * type_nice) +
              (Heading("Diverse") * diverse)*(Heading("Actor type") * type_nice) +
              (Heading("University education") * university_ed)*(Heading("Actor type") * type_nice) +
              (Heading("Employed") * employed1)*(Heading("Actor type") * type_nice) +
              (Heading("Workers treatment") * treatment_workers)*(Heading("Actor type") * type_nice) ~
              Mean + SD + Min + Median + Max + N,
            data=df3, escape = FALSE, output="tables/TableA1.tex",
            title="Summary statistics by actor type ⁠\\label{tab:descriptive}⁠")

vtable::st(df3)} else{
  
  message("To protect respondents privacy, the age data of legislators is excluded. It is removed from this table.")
  
  datasummary((Heading("Infrastructure projects") * infrastructure)*(Heading("Actor type") * type_nice) +
                (Heading("Subsidies") * subsidies)*(Heading("Actor type") * type_nice) +
                (Heading("Tariff increases") * tariffs)*(Heading("Actor type") * type_nice) +
                (Heading("Tax cuts") * tax)*(Heading("Actor type") * type_nice) +
                #Heading("Legislator") * legislator)*(Heading("Actor type") * type) +
                (Heading("Political ideology") * ideology)*(Heading("Actor type") * type_nice) +
                (Heading("Male") * male)*(Heading("Actor type") * type_nice) +
                (Heading("Diverse") * diverse)*(Heading("Actor type") * type_nice) +
                (Heading("University education") * university_ed)*(Heading("Actor type") * type_nice) +
                (Heading("Employed") * employed1)*(Heading("Actor type") * type_nice) +
                (Heading("Workers treatment") * treatment_workers)*(Heading("Actor type") * type_nice) ~
                Mean + SD + Min + Median + Max + N,
              data=df3, escape = FALSE, output="tables/TableA1.tex",
              title="Summary statistics by actor type ⁠\\label{tab:descriptive}⁠")
}

# Descriptive regression ####

df3$subsidies <- factor(df3$subsidies)
df3$tariffs <- factor(df3$tariffs)
df3$tax <- factor(df3$tax)
df3$infrastructure <- factor(df3$infrastructure)
df3$type <- factor(df3$type)
df3$country <- factor(df3$country)

# .. Table B1 ####


if (dir.exists("not_sharseable/")) {

m1sl <- clmm2(location = subsidies ~ age + gender + ideology + university + treatment, random = country, Hess = T, data= df3, subset = type == "legislator")
m1tl <- clmm2(location = tariffs ~ age + gender + ideology + university + treatment, random = country, Hess = T, data= df3, subset = type == "legislator")
m1xl <- clmm2(location = tax ~ age + gender + ideology + university + treatment, random = country, Hess = T, data= df3, subset = type == "legislator")
m1il <- clmm2(location = infrastructure ~ age + gender + ideology + university + treatment, random = country, Hess = T, data= df3, subset = type == "legislator")

modelsummary::modelsummary(models = list("Infrastructure" = m1il, "Subsidies" = m1sl, "Tariffs" = m1tl, "Tax cuts" = m1xl),
                           output = "tables/TableB1.tex",
                           title = "Legislators and predictors of policy support",
                           label = "tab:model1_leg",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('age' = 'Age',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'ideology' = 'Ideology',
                                        'universityyes' = 'University (ref. no university)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))} else{
  
  message("To protect respondents privacy, the age data of legislators is excluded. It is removed from this table.")
                           }

# .. Table B2 ####


if (dir.exists("not_shareable/")) {
  m1sc <- clmm2(location = subsidies ~ age + gender + ideology + university + employ + treatment, random = country, Hess = T, data= df3, subset = type == "public")
  m1tc <- clmm2(location = tariffs ~ age + gender + ideology + university + employ + treatment, random = country, Hess = T, data= df3, subset = type == "public")
  m1xc <- clmm2(location = tax ~ age + gender + ideology + university + employ + treatment, random = country, Hess = T, data= df3, subset = type == "public")
  m1ic <- clmm2(location = infrastructure ~ age + gender + ideology + university + employ + treatment, random = country, Hess = T, data= df3, subset = type == "public")
  
  modelsummary::modelsummary(models = list("Infrastructure" = m1ic, "Subsidies" = m1sc, "Tariffs" = m1tc, "Tax cuts" = m1xc),
                             output = "tables/TableB2.tex",
                             title = "The public and predictors of policy support",
                             label = "tab:model1_pub",
                             estimate = "{estimate}{stars}", 
                             stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                             notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                             fmt = 2,
                             coef_map = c('age' = 'Age',
                                          'genderMale' = 'Male (ref. female)',
                                          'genderOther' = 'Other (ref. female)',
                                          'ideology' = 'Ideology',
                                          'universityyes' = 'University (ref. no university)',
                                          'employEmployed' = 'Employed (ref. unemployed)',
                                          'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                          '0|1' = '0|1',
                                          '1|2' = '1|2',
                                          '2|3' = '2|3',
                                          '3|4' = '3|4'),
                             gof_map = c("nobs", "aic", "bic"))} else{
                               
                               message("To protect respondents privacy, the age data of legislators is excluded. It is removed from this table.")
                             }

## .. Table B3 ####

if (dir.exists("not_shareable/")) {
  
  m1slc <- clmm2(location = subsidies ~ (age + gender + ideology + university + treatment) * type, random = country, Hess = T, data= df3)
  m1tlc <- clmm2(location = tariffs ~ (age + gender + ideology + university + treatment) * type, random = country, Hess = T, data= df3)
  m1xlc <- clmm2(location = tax ~ (age + gender + ideology + university + treatment) * type, random = country, Hess = T, data= df3)
  m1ilc <- clmm2(location = infrastructure ~ (age + gender + ideology + university + treatment) * type, random = country, Hess = T, data= df3)
  
  modelsummary::modelsummary(models = list("Infrastructure" = m1ilc, "Subsidies" = m1slc, "Tariffs" = m1tlc, "Tax cuts" = m1xlc),
                             output = "tables/TableB3.tex",
                             title = "Differences in predictors of policy support",
                             label = "tab:model1_interaction",
                             estimate = "{estimate}{stars}", 
                             stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                             notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                             fmt = 2,
                             coef_map = c('age' = 'Age',
                                          'genderMale' = 'Male (ref. female)',
                                          'genderOther' = 'Other (ref. female)',
                                          'ideology' = 'Ideology',
                                          'universityyes' = 'University (ref. no university)',
                                          'employEmployed' = 'Employed (ref. unemployed)',
                                          'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                          'typepublic' = 'Public (ref. legislators)',
                                          'age:typepublic' = "Age × Public",
                                          'genderMale:typepublic' = "Male × Public",
                                          'ideology:typepublic' = "Ideology × Public",
                                          'universityyes:typepublic' = "University × Public",
                                          'treatmentworkers:typepublic' = "Worker × Public",
                                          '0|1' = '0|1',
                                          '1|2' = '1|2',
                                          '2|3' = '2|3',
                                          '3|4' = '3|4'),
                             gof_map = c("nobs", "aic", "bic"))} else{
                               
                               message("To protect respondents privacy, the age data of legislators is excluded. It is removed from this table.")
                             }


# .. Table B4 ####

df3$type3 <- factor(ifelse(df3$type == "public" & df3$university == "yes", "University public",
                           ifelse(df3$type == "public" & df3$university == "no", "Non-university public",
                                  ifelse(df3$type == "legislator", "Legislator", NA))),
                    levels = c("Legislator", "University public", "Non-university public"))


m1s3 <- clmm2(location = subsidies ~ type3 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1t3 <- clmm2(location = tariffs ~ type3 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1x3 <- clmm2(location = tax ~ type3 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1i3 <- clmm2(location = infrastructure ~ type3 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)

modelsummary::modelsummary(models = list("Infrastructure" = m1i3, "Subsidies" = m1s3, "Tariffs" = m1t3, "Tax cuts" = m1x3),
                           output = "tables/TableB4.tex",
                           title = "Legislators, the public, education and policy support",
                           label = "tab:model1education",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('type3Non-university public' = 'Non-univ. educated Public (ref. legislators)',
                                        'type3University public' = 'Univ. educated Public (ref. legislators)',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'ideology' = 'Ideology',
                                        'universityyes' = 'University (no university)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))


# .. Table B5 ####

df3$type2 <- factor(ifelse(df3$type == "public" & df3$employ == "Employed", "Employed public",
                           ifelse(df3$type == "public" & df3$employ == "Unemployed", "Unemployed public",
                                  ifelse(df3$type == "legislator", "Legislator", NA))),
                    levels = c("Legislator", "Employed public", "Unemployed public"))

m1s2 <- clmm2(location = subsidies ~ type2 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1t2 <- clmm2(location = tariffs ~ type2 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1x2 <- clmm2(location = tax ~ type2 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1i2 <- clmm2(location = infrastructure ~ type2 + gender + ideology + university + treatment, random = country, Hess = T, data= df3)

modelsummary::modelsummary(models = list("Infrastructure" = m1i2, "Subsidies" = m1s2, "Tariffs" = m1t2, "Tax cuts" = m1x2),
                           output = "tables/TableB5.tex",
                           title = "Legislators, the public, employment status and policy support",
                           label = "tab:model1employment",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('type2Unemployed public' = 'Unemployed Public (ref. legislators)',
                                        'type2Employed public' = 'Employed Public (ref. legislators)',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'ideology' = 'Ideology',
                                        'universityyes' = 'University (no university)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# Main analyses ####


m1s <- clmm2(location = subsidies ~ type + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1t <- clmm2(location = tariffs ~ type + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1x <- clmm2(location = tax ~ type + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m1i <- clmm2(location = infrastructure ~ type + gender + ideology + university + treatment, random = country, Hess = T, data= df3)

# .. Calculation of effect sizes for discussion in paper ####

marginaleffects::predictions(m1s, newdata = "mean", by = c("subsidies", "type")) %>% 
  data.frame() %>% 
  select(subsidies, type, estimate) %>% 
  arrange(subsidies) %>% 
  group_by(subsidies) %>% 
  summarise(leg_pub = estimate[1]-estimate[2])

marginaleffects::predictions(m1t, newdata = "mean", by = c("tariffs", "type")) %>% 
  data.frame() %>% 
  select(tariffs, type, estimate) %>% 
  arrange(tariffs) %>% 
  group_by(tariffs) %>% 
  summarise(leg_pub = estimate[1]-estimate[2])

marginaleffects::avg_predictions(m1x, newdata = "mean", by = c("tax", "type")) %>% 
  data.frame() %>% 
  select(tax, type, estimate) %>% 
  arrange(tax) %>% 
  group_by(tax) %>% 
  summarise(leg_pub = estimate[1]-estimate[2])

marginaleffects::avg_predictions(m1i, newdata = "mean", by = c("infrastructure", "type")) %>% 
  data.frame() %>% 
  select(infrastructure, type, estimate) %>% 
  arrange(infrastructure) %>% 
  group_by(infrastructure) %>% 
  summarise(leg_pub = estimate[1]-estimate[2])

# .. Table C1 ####

modelsummary::modelsummary(models = list("Infrastructure" = m1i, "Subsidies" = m1s, "Tariffs" = m1t, "Tax cuts" = m1x),
                           output = "tables/TableC1.tex",
                           title = "Legislators, the public and policy support",
                           label = "tab:model1",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('typepublic' = 'Public (ref. legislators)',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'ideology' = 'Ideology',
                                        'universityyes' = 'University (no university)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# .. Figure 4 ####


data.frame(out = factor(c("Subsidies", "Tariffs", "Taxes", "Infrastructure"),
                        levels = c("Infrastructure", "Subsidies", "Tariffs", "Taxes"),
                        labels = c("Infrastructure projects", "Subsidies", "Tariffs", "Tax cuts")),
           coef = c(m1s$coefficients["typepublic"],
                    m1t$coefficients["typepublic"],
                    m1x$coefficients["typepublic"],
                    m1i$coefficients["typepublic"]),
           se = c(summary(m1s)$coefficients["typepublic",2],
                  summary(m1t)$coefficients["typepublic",2],
                  summary(m1x)$coefficients["typepublic",2],
                  summary(m1i)$coefficients["typepublic",2])) %>% 
  ggplot(., aes(x = out, y = coef, ymin = coef - 1.96*se, ymax = coef + 1.96*se)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  geom_linerange(aes(ymin = coef - 1.645*se, ymax = coef + 1.645*se), linewidth = 0.75, position = position_dodge(width = 1)) +
  theme_minimal() +
  labs(x = "Compensatory measure", y = "Coefficient for public (vs legislators)") + 
  theme(legend.position = "bottom") + 
  NULL

ggsave("figures/Figure4.pdf",
       width = 16, height = 9, units = "cm")

# Interaction with ideology ####

m2s <- clmm2(location = subsidies ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m2t <- clmm2(location = tariffs ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m2x <- clmm2(location = tax ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data= df3)
m2i <- clmm2(location = infrastructure ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data= df3)

# .. Table C2 ####

modelsummary::modelsummary(models = list("Infrastructure" = m2i, "Subsidies" = m2s, "Tariffs" = m2t, "Tax cuts" = m2x),
                           output = "tables/TableC2.tex",
                           title = "Ideology, actor type and policy support",
                           label = "tab:model2",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('typepublic' = 'Public (ref. legislators)',
                                        'ideology' = 'Ideology',
                                        'typepublic:ideology' = 'Public \\times Ideology',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'universityyes' = 'University (no university)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# .. Figure C1 ####

#If the bootstrapping data already exists, the loop is jumped

if (!file.exists("output/bootstrapping_margins_m2.rds")) {
  
  out_df_me <- data.frame(out = as.character(),
                          dv = as.character(),
                          term = as.character(),
                          contrast = as.character(),
                          ideology = as.numeric(),
                          estimate = as.numeric(),
                          id = as.numeric())
  
  #select(out, dv, term, contrast, ideology, estimate),
  
  nboot <- 1000
  run_start <- Sys.time()
  set.seed(123)
  
  pb = txtProgressBar(min = 0, max = nboot, initial = 0, style = 3) 
  
  for (i in 1:nboot){
    
    temp_df <- df3[sample(nrow(df3), nrow(df3), replace = T), ]
    
    m2st <- clmm2(location = subsidies ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2tt <- clmm2(location = tariffs ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2xt <- clmm2(location = tax ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2it <- clmm2(location = infrastructure ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    
    out_temp <- rbind(
      marginaleffects::slopes(m2st, variable = "type", by = c("subsidies", "ideology")) %>% 
        mutate(out = "Subsidies") %>% 
        rename(dv = subsidies) %>% 
        select(out, dv, term, contrast, ideology, estimate),
      marginaleffects::slopes(m2tt, variable = "type", by = c("tariffs", "ideology")) %>% 
        mutate(out = "Tariffs") %>% 
        rename(dv = tariffs) %>% 
        select(out, dv, term, contrast, ideology, estimate),
      marginaleffects::slopes(m2xt, variable = "type", by = c("tax", "ideology")) %>%
        mutate(out = "Tax cuts") %>% 
        rename(dv = tax) %>% 
        select(out, dv, term, contrast, ideology, estimate),
      marginaleffects::slopes(m2it, variable = "type", by = c("infrastructure", "ideology")) %>%         mutate(out = "Infrastructure") %>% 
        rename(dv = infrastructure) %>% 
        select(out, dv, term, contrast, ideology, estimate)) %>% 
      mutate(id = i,
             dv = as.character(dv))
    
    out_df_me <- bind_rows(out_df_me,
                           out_temp)
    
    setTxtProgressBar(pb,i)
    
    
  }
  
  run_end <- Sys.time()
  run_end - run_start
  
  write_rds(out_df_me, "output/bootstrapping_margins_m2.rds")
} else{
  out_df_me <- readRDS("output/bootstrapping_margins_m2.rds")
}

out_df_me %>% 
  group_by(out, dv, ideology) %>% 
  summarise(lower = quantile(estimate, probs = 0.025, na.rm=T),
            upper = quantile(estimate, probs = 0.975, na.rm=T),
            estimate = mean(estimate, na.rm=T),
            sign = ifelse(upper < 0 | lower > 0, "*", "n.s.")) %>% 
  mutate(dv = factor(dv,
                     levels = c(0:4),
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))) %>% 
  ggplot(., aes(x = ideology, y = estimate, ymin = lower, ymax = upper, lty = sign)) +
  geom_hline(yintercept = 0, lty = "dotted") +
  geom_point() +
  geom_errorbar(width = 0.1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_grid(out~dv, scales = "free_y") +
  labs(x = "Ideology", y = "Coefficient for public (vs legislators)") +
  scale_linetype_discrete("Significance") +
  scale_y_continuous(minor_breaks = c(-.3,-.2,-.1, 0, .1, .2, 0.3)) +
  scale_x_continuous(minor_breaks = c(1:5)) +
  NULL

ggsave("figures/FigureC1.pdf",
       width = 21, height = 16, units = "cm")

# .. Figure 5 ####

if (!file.exists("output/bootstrapping_predictions_model2.rds")) {
  
  out_df_pred <- data.frame(out = as.character(),
                            dv = as.character(),
                            type = as.character(),
                            ideology = as.numeric(),
                            estimate = as.numeric(), 
                            id = as.numeric())
  
  #select(out, dv, term, contrast, ideology, estimate),
  
  nboot <- 1000
  run_start <- Sys.time()
  set.seed(123)
  
  pb = txtProgressBar(min = 0, max = nboot, initial = 0, style = 3) 
  
  for (i in 1:nboot){
    
    temp_df <- df3[sample(nrow(df3), nrow(df3), replace = T), ]
    
    m2st <- clmm2(location = subsidies ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2tt <- clmm2(location = tariffs ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2xt <- clmm2(location = tax ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2it <- clmm2(location = infrastructure ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    
    out_temp <- rbind(
      marginaleffects::predictions(m2st, by = c("subsidies", "ideology", "type")) %>% 
        mutate(out = "Subsidies") %>% 
        select(out, dv = subsidies, type, ideology, estimate),
      marginaleffects::predictions(m2tt, by = c("tariffs", "ideology", "type")) %>% 
        mutate(out = "Tariffs") %>% 
        select(out, dv = tariffs, type, ideology, estimate),
      marginaleffects::predictions(m2xt, by = c("tax", "ideology", "type")) %>% 
        mutate(out = "Tax cuts") %>% 
        select(out, dv = tax, type, ideology, estimate),
      marginaleffects::predictions(m2it, by = c("infrastructure", "ideology", "type")) %>% 
        mutate(out = "Infrastructure") %>% 
        select(out, dv = infrastructure , type, ideology, estimate)) %>% 
      mutate(id = i)
    
    out_df_pred <- bind_rows(out_df_pred,
                             out_temp)
    
    setTxtProgressBar(pb,i)
    
  }
  
  run_end <- Sys.time()
  run_end - run_start
  
  write_rds(out_df_pred, "output/bootstrapping_predictions_model2.rds")
} else{
  out_df_pred <- readRDS("output/bootstrapping_predictions_model2.rds")
}


out_df_pred %>% 
  group_by(out, dv, ideology, type) %>% 
  summarise(lower = quantile(estimate, probs = 0.025, na.rm=T),
            upper = quantile(estimate, probs = 0.975, na.rm=T),
            estimate = mean(estimate, na.rm=T)) %>% 
  mutate(dv = factor(dv,
                     levels = c(0:4),
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))) %>% 
  ggplot(., aes(x = ideology, y = estimate, ymin = lower, ymax = upper, fill = type)) +
  geom_ribbon(alpha = .3) +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_grid(out~dv, scales = "free_y") +
  labs(x = "Ideology", y = "Predicted support") +
  scale_fill_discrete("Type", labels = c("Legislators", "Public")) +
  scale_y_continuous(minor_breaks = c(0, 0.2, 0.4, 0.6)) +
  scale_x_continuous(minor_breaks = c(1:5)) +
  NULL

ggsave("figures/Figure5.pdf",
       width = 21, height = 16, units = "cm")

# Subgroup analyses for the UK and Spain ####

df_ukes <- df3 %>%  filter(country %in% c("Spain", "United Kingdom")) %>% mutate(country = droplevels(country, except = c("Spain", "United Kingdom")))

# .. Table E1 ####


m1ssub <- clm(subsidies ~ type + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)
m1tsub <- clm(tariffs ~ type + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)
m1xsub <- clm(tax ~ type + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)
m1isub <- clm(infrastructure ~ type + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)

modelsummary::modelsummary(models = list("Infrastructure" = m1isub, "Subsidies" = m1ssub, "Tariffs" = m1tsub, "Tax cuts" = m1xsub),
                           output = "tables/TableE1.tex",
                           title = "Legislators, the public and policy support",
                           label = "tab:model1_UKES",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('typepublic' = 'Public (ref. legislators)',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'ideology' = 'Ideology',
                                        'universityyes' = 'University (no university)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        'countryUnited Kingdom' = 'UK (ref. ES)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# .. Table E2 ####


m2ssub <- clm(formula = subsidies ~ type * ideology + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)
m2tsub <- clm(formula = tariffs ~ type * ideology + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)
m2xsub <- clm(formula = tax ~ type * ideology + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)
m2isub <- clm(formula = infrastructure ~ type * ideology + gender + ideology + university + treatment + country, Hess = T, data = df_ukes)


modelsummary::modelsummary(models = list("Infrastructure" = m2isub, "Subsidies" = m2ssub, "Tariffs" = m2tsub, "Tax cuts" = m2xsub),
                           output = "tables/TableE2.tex",
                           title = "Ideology, actor type and policy support",
                           label = "tab:model2_UKES",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('typepublic' = 'Public (ref. legislators)',
                                        'ideology' = 'Ideology',
                                        'typepublic:ideology' = 'Public \times Ideology',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'universityyes' = 'University (no university)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        'countryUnited Kingdom' = 'UK (ref. ES)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# .. Figure 6 ####

data.frame(out = factor(rep(c("Subsidies", "Tariffs", "Taxes", "Infrastructure"),2),
                        levels = c("Infrastructure", "Subsidies", "Tariffs", "Taxes"),
                        labels = c("Infrastructure projects", "Subsidies", "Tariffs", "Tax cuts")),
           coef = c(m1s$coefficients["typepublic"],
                    m1t$coefficients["typepublic"],
                    m1x$coefficients["typepublic"],
                    m1i$coefficients["typepublic"],
                    m1ssub$coefficients["typepublic"],
                    m1tsub$coefficients["typepublic"],
                    m1xsub$coefficients["typepublic"],
                    m1isub$coefficients["typepublic"]),
           se = c(summary(m1s)$coefficients["typepublic",2],
                  summary(m1t)$coefficients["typepublic",2],
                  summary(m1x)$coefficients["typepublic",2],
                  summary(m1i)$coefficients["typepublic",2],
                  summary(m1ssub)$coefficients["typepublic",2],
                  summary(m1tsub)$coefficients["typepublic",2],
                  summary(m1xsub)$coefficients["typepublic",2],
                  summary(m1isub)$coefficients["typepublic",2]),
           model = factor(rep(c("Main models", "UK/Spain subset"), each = 4),
                          levels = c("Main models", "UK/Spain subset"))) %>% 
  ggplot(., aes(x = out, y = coef,
                ymin = coef - 1.96*se, ymax = coef + 1.96*se,
                colour = model, shape = model)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_linerange(aes(ymin = coef - 1.645*se, ymax = coef + 1.645*se), linewidth = 0.75, position = position_dodge(width = .5)) +
  theme_minimal() +
  labs(x = "Compensatory measure", y = "Coefficient for public (vs legislators)") + 
  scale_colour_discrete("Model") +
  scale_shape_discrete("Model") +
  theme(legend.position = "bottom") + 
  NULL

ggsave("figures/Figure6.pdf",
       width = 16, height = 9, units = "cm")

# Jackknifing ####

#If jackknifing data already exists, the loop is jumped

if (!file.exists("output/jacknifing.rds")) {
  
  country_list <- unique(df3$country)
  
  df_out <- data.frame(country = as.character(),
                       dv = as.character(),
                       mod = as.double(),
                       var = as.character(),
                       est = as.numeric(),
                       se = as.numeric())
  
  j <- length(country_list)
  
  for (i in country_list){
    
    temp_dat <- df3 %>% filter(country != i)
    
    m1stemp <- clmm2(location = subsidies ~ type + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    m1ttemp <- clmm2(location = tariffs ~ type + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    m1xtemp <- clmm2(location = tax ~ type + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    m1itemp <- clmm2(location = infrastructure ~ type + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    
    
    m2stemp <- clmm2(location = subsidies ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    m2ttemp <- clmm2(location = tariffs ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    m2xtemp <- clmm2(location = tax ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    m2itemp <- clmm2(location = infrastructure ~ type * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_dat)
    
    temp_out <- data.frame(country = rep(i, 12),
                           dv = c("Subsidies", "Tariffs", "Tax cut","Infrastructure", rep(c("Subsidies", "Tariffs", "Tax cut","Infrastructure"),each = 2)),
                           mod = c(rep(c(1,2,2), each = 4)),
                           var = c(rep("type", 4), rep(c("type", "type:ideology"), 4)),
                           est = c(summary(m1stemp)$coefficients["typepublic",1],
                                   summary(m1ttemp)$coefficients["typepublic",1],
                                   summary(m1xtemp)$coefficients["typepublic",1],
                                   summary(m1itemp)$coefficients["typepublic",1],
                                   summary(m2stemp)$coefficients["typepublic",1],
                                   summary(m2stemp)$coefficients["typepublic:ideology",1],
                                   summary(m2ttemp)$coefficients["typepublic",1],
                                   summary(m2ttemp)$coefficients["typepublic:ideology",1],
                                   summary(m2xtemp)$coefficients["typepublic",1],
                                   summary(m2xtemp)$coefficients["typepublic:ideology",1],
                                   summary(m2itemp)$coefficients["typepublic",1],
                                   summary(m2itemp)$coefficients["typepublic:ideology",1]),
                           se = c(summary(m1stemp)$coefficients["typepublic",2],
                                  summary(m1ttemp)$coefficients["typepublic",2],
                                  summary(m1xtemp)$coefficients["typepublic",2],
                                  summary(m1itemp)$coefficients["typepublic",2],
                                  summary(m2stemp)$coefficients["typepublic",2],
                                  summary(m2stemp)$coefficients["typepublic:ideology",2],
                                  summary(m2ttemp)$coefficients["typepublic",2],
                                  summary(m2ttemp)$coefficients["typepublic:ideology",2],
                                  summary(m2xtemp)$coefficients["typepublic",2],
                                  summary(m2xtemp)$coefficients["typepublic:ideology",2],
                                  summary(m2itemp)$coefficients["typepublic",2],
                                  summary(m2itemp)$coefficients["typepublic:ideology",2]))
    
    df_out <- bind_rows(df_out,
                        temp_out)
    
    j <- j - 1
    
    print(paste0(i, ": Done! ", j, " to go."))
    
  }
  
  saveRDS(df_out, file = "output/jacknifing.rds")} else{
    df_out <- readRDS("output/jacknifing.rds")
  }

# .. Figure F1 ####

df_out %>% 
  filter(country != "United Kingdom" & country != "Spain" & country != "Poland") %>% 
  mutate(var = factor(ifelse(var == "type", "Public (ref. legislators)", "Interaction term"),
                      levels = c("Public (ref. legislators)", "Interaction term")),
         country_short = countrycode::countrycode(country, "country.name", "iso3c"),
         qoi = paste0("Model ", mod, " - ", var)) %>% 
  ggplot(., aes(x = forcats::fct_rev(country_short), y = est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_pointrange(size = .25) +
  facet_grid(dv ~ qoi) +
  coord_flip() +
  geom_hline(yintercept = 0, lty = "dotted") +
  theme_minimal() + 
  labs(x = "Country", y = "Coefficient") +
  NULL

ggsave("figures/FigureF1.pdf",
       width = 21, height = 25, units = "cm")

# Treatment condition ####

# .. Table G1 ####

m1streat <- clmm2(subsidies ~ type * treatment + gender  + ideology + university + treatment, random = country, Hess = T, data = df3)
m1ttreat <- clmm2(tariffs ~ type * treatment + gender  + ideology + university + treatment, random = country, Hess = T, data = df3)
m1xtreat <- clmm2(tax ~ type * treatment + gender  + ideology + university + treatment, random = country, Hess = T, data = df3)
m1itreat <- clmm2(infrastructure ~ type * treatment + gender  + ideology + university + treatment, random = country, Hess = T, data = df3)


modelsummary::modelsummary(models = list("Infrastructure" = m1itreat, "Subsidies" = m1streat, "Tariffs" = m1ttreat, "Tax cuts" = m1xtreat),
                           output = "tables/TableG1.tex",
                           title = "Legislators, the public, treatment and policy support",
                           label = "tab:model1treatment",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('typepublic' = 'Public (ref. legislators)',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        'typepublic:treatmentworkers' = 'Public \\times Worker',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'ideology' = 'Ideology',
                                        'universityyes' = 'University (no university)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# .. Table G2 ####


m2streat <- clmm2(subsidies ~ type * ideology * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = df3)
m2ttreat <- clmm2(tariffs ~ type * ideology * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = df3)
m2xtreat <- clmm2(tax ~ type * ideology * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = df3)
m2itreat <- clmm2(infrastructure ~ type * ideology * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = df3)

modelsummary::modelsummary(models = list("Infrastructure" = m2itreat, "Subsidies" = m2streat, "Tariffs" = m2ttreat, "Tax cuts" = m2xtreat),
                           output = "tables/TableG2.tex",
                           title = "Legislators, the public, treatment and policy support",
                           label = "tab:model1treatment",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('typepublic' = 'Public (ref. legislators)',
                                        'ideology' = 'Ideology',
                                        'treatmentworkers' = 'Worker Treatment (ref. firms)',
                                        'typepublic:ideology' = 'Public \\times Ideology',
                                        'typepublic:treatmentworkers' = 'Public \\times Worker',
                                        'ideology:treatmentworkers' = 'Ideology \\times Worker',
                                        'typepublic:ideology:treatmentworkers' = 'Public \\times Ideology \\times Worker',
                                        'genderMale' = 'Male (ref. female)',
                                        'genderOther' = 'Other (ref. female)',
                                        'universityyes' = 'University (no university)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# .. Figure G1 ####


if (!file.exists("output/bootstrapping_margins_m1_treatment.rds")) {
  
  out_df_me_treatment <- data.frame(out = as.character(),
                                    dv = as.character(),
                                    term = as.character(),
                                    contrast = as.character(),
                                    type = as.numeric(),
                                    estimate = as.numeric(),
                                    id = as.numeric())
  nboot <- 1000
  run_start <- Sys.time()
  set.seed(123)
  
  pb = txtProgressBar(min = 0, max = nboot, initial = 0, style = 3) 
  
  for (i in 1:nboot){
    
    temp_df <- df3[sample(nrow(df3), nrow(df3), replace = T), ]
    
    m1st <- clmm2(location = subsidies ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m1tt <- clmm2(location = tariffs ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m1xt <- clmm2(location = tax ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m1it <- clmm2(location = infrastructure ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    
    out_temp <- rbind(
      marginaleffects::slopes(m1st, variable = "type", by = c("subsidies", "treatment")) %>% 
        mutate(out = "Subsidies") %>% 
        rename(dv = subsidies) %>% 
        select(out, dv, term, contrast, treatment, estimate),
      marginaleffects::slopes(m1tt, variable = "type", by = c("tariffs", "treatment")) %>% 
        mutate(out = "Tariffs") %>% 
        rename(dv = tariffs) %>% 
        select(out, dv, term, contrast, treatment, estimate),
      marginaleffects::slopes(m1xt, variable = "type", by = c("tax", "treatment")) %>%
        mutate(out = "Tax cuts") %>% 
        rename(dv = tax) %>% 
        select(out, dv, term, contrast, treatment, estimate),
      marginaleffects::slopes(m1it, variable = "type", by = c("infrastructure", "treatment")) %>% 
        mutate(out = "Infrastructure") %>% 
        rename(dv = infrastructure) %>% 
        select(out, dv, term, contrast, treatment, estimate)) %>% 
      mutate(id = i,
             dv = as.character(dv))
    
    out_df_me_treatment <- bind_rows(out_df_me_treatment,
                                     out_temp)
    
    setTxtProgressBar(pb,i)
    
    
  }
  
  run_end <- Sys.time()
  run_end - run_start
  
  write_rds(out_df_me_treatment, "output/bootstrapping_margins_m1_treatment.rds")
} else{
  out_df_me_treatment <- readRDS("output/bootstrapping_margins_m1_treatment.rds")
}


out_df_me_treatment %>% 
  group_by(out, dv, treatment) %>% 
  summarise(lower = quantile(estimate, probs = 0.025, na.rm=T),
            upper = quantile(estimate, probs = 0.975, na.rm=T),
            estimate = mean(estimate, na.rm=T),
            sign = ifelse(upper < 0 | lower > 0, "*", "n.s.")) %>% 
  mutate(dv = factor(dv,
                     levels = c(0:4),
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))) %>% 
  ggplot(., aes(x = treatment, y = estimate, ymin = lower, ymax = upper, lty = sign)) +
  geom_hline(yintercept = 0, lty = "dotted") +
  geom_point(position = position_dodge(width = 1), ) +
  geom_errorbar(position = position_dodge(width = 1), width = 0.1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_grid(out~dv, scales = "free_y") +
  labs(x = "Ideology", y = "Marginal effect of public (vs legislators)") +
  scale_linetype_discrete("Significance") +
  scale_colour_discrete("Treatment", labels = c("Firms", "Worker")) + 
  scale_shape_discrete("Treatment", labels = c("Firms", "Worker")) + 
  NULL

ggsave("figures/FigureG1.pdf",
       width = 21, height = 16, units = "cm")

# .. Figure G2 ####


if (!file.exists("output/bootstrapping_predictions_m1_treatment.rds")) {
  
  out_df_me_pred_treatment <- data.frame(out = as.character(),
                                         dv = as.character(),
                                         type = as.character(),
                                         treatment = as.character(),
                                         estimate = as.numeric(),
                                         id = as.numeric())
  
  nboot <- 1000
  run_start <- Sys.time()
  set.seed(123)
  
  pb = txtProgressBar(min = 0, max = nboot, initial = 0, style = 3) 
  
  for (i in 1:nboot){
    
    temp_df <- df3[sample(nrow(df3), nrow(df3), replace = T), ]
    
    m1st <- clmm2(location = subsidies ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m1tt <- clmm2(location = tariffs ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m1xt <- clmm2(location = tax ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m1it <- clmm2(location = infrastructure ~ type * treatment + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    
    out_temp <- rbind(
      marginaleffects::predictions(m1st, by = c("subsidies", "type", "treatment")) %>% 
        mutate(out = "Subsidies") %>% 
        rename(dv = subsidies) %>% 
        select(out, dv, type, treatment, estimate),
      marginaleffects::predictions(m1tt, by = c("tariffs", "type", "treatment")) %>% 
        mutate(out = "Tariffs") %>% 
        rename(dv = tariffs) %>% 
        select(out, dv, type, treatment, estimate),
      marginaleffects::predictions(m1xt, by = c("tax", "type", "treatment")) %>% 
        mutate(out = "Tax cuts") %>% 
        rename(dv = tax) %>% 
        select(out, dv, type, treatment, estimate),
      marginaleffects::predictions(m1it, by = c("infrastructure", "type", "treatment")) %>% 
        mutate(out = "Infrastructure") %>% 
        rename(dv = infrastructure) %>% 
        select(out, dv, type, treatment, estimate)) %>% 
      mutate(id = i,
             dv = as.character(dv))
    
    out_df_me_pred_treatment <- bind_rows(out_df_me_pred_treatment,
                                          out_temp)
    
    setTxtProgressBar(pb,i)
    
    
  }
  
  run_end <- Sys.time()
  run_end - run_start
  
  write_rds(out_df_me_pred_treatment, "output/bootstrapping_predictions_m1_treatment.rds")
} else{
  out_df_me_pred_treatment <- readRDS("output/bootstrapping_predictions_m1_treatment.rds")
}


out_df_me_pred_treatment %>% 
  group_by(out, dv, type, treatment) %>% 
  summarise(lower = quantile(estimate, probs = 0.025, na.rm=T),
            upper = quantile(estimate, probs = 0.975, na.rm=T),
            estimate = mean(estimate, na.rm=T)) %>% 
  mutate(dv = factor(dv,
                     levels = c(0:4),
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))) %>% 
  ggplot(., aes(x = treatment, y = estimate, ymin = lower, ymax = upper, colour = type, shape = type)) +
  geom_point(position = position_dodge(width = 1), ) +
  geom_errorbar(position = position_dodge(width = 1), width = 0.1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_grid(out~dv, scales = "free_y") +
  labs(x = "Treatment", y = "Predicted probability") +
  scale_y_continuous(minor_breaks = c(0, .1, .2, 0.3, 0.4, 0.5)) +
  scale_x_discrete(labels = c("Firms", "Worker")) + 
  scale_colour_discrete("Type", labels = c("Legislator", "Public")) + 
  scale_shape_discrete("Type", labels = c("Legislator", "Public")) + 
  NULL

ggsave("figures/FigureG2.pdf",
       width = 21, height = 16, units = "cm")

# .. Figure G3 ####


if (!file.exists("output/bootstrapping_margins_m2_treatment.rds")) {
  
  out_df_me_treatmentideology <- data.frame(out = as.character(),
                                            dv = as.character(),
                                            term = as.character(),
                                            contrast = as.character(),
                                            ideology = as.numeric(),
                                            estimate = as.numeric(),
                                            id = as.numeric())
  
  nboot <- 1000
  run_start <- Sys.time()
  set.seed(123)
  
  pb = txtProgressBar(min = 0, max = nboot, initial = 0, style = 3) 
  
  for (i in 1:nboot){
    
    temp_df <- df3[sample(nrow(df3), nrow(df3), replace = T), ]
    
    m2st <- clmm2(location = subsidies ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2tt <- clmm2(location = tariffs ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2xt <- clmm2(location = tax ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2it <- clmm2(location = infrastructure ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    
    out_temp <- rbind(
      marginaleffects::slopes(m2st, variable = "type", by = c("subsidies", "ideology", "treatment")) %>% 
        mutate(out = "Subsidies") %>% 
        rename(dv = subsidies) %>% 
        select(out, dv, term, contrast, treatment, ideology, estimate),
      marginaleffects::slopes(m2tt, variable = "type", by = c("tariffs", "ideology", "treatment")) %>% 
        mutate(out = "Tariffs") %>% 
        rename(dv = tariffs) %>% 
        select(out, dv, term, contrast, treatment, ideology, estimate),
      marginaleffects::slopes(m2xt, variable = "type", by = c("tax", "ideology", "treatment")) %>%
        mutate(out = "Tax cuts") %>% 
        rename(dv = tax) %>% 
        select(out, dv, term, contrast, treatment, ideology, estimate),
      marginaleffects::slopes(m2it, variable = "type", by = c("infrastructure", "ideology", "treatment")) %>% 
        mutate(out = "Infrastructure") %>% 
        rename(dv = infrastructure) %>% 
        select(out, dv, term, contrast, treatment, ideology, estimate)) %>% 
      mutate(id = i,
             dv = as.character(dv))
    
    out_df_me_treatmentideology <- bind_rows(out_df_me_treatmentideology,
                                             out_temp)
    
    setTxtProgressBar(pb,i)
    
    
  }
  
  run_end <- Sys.time()
  run_end - run_start
  
  write_rds(out_df_me_treatmentideology, "output/bootstrapping_margins_m2_treatment.rds")
} else{
  out_df_me_treatmentideology <- readRDS("output/bootstrapping_margins_m2_treatment.rds")
}


out_df_me_treatmentideology %>% 
  group_by(out, dv, ideology, treatment) %>% 
  summarise(lower = quantile(estimate, probs = 0.025, na.rm=T),
            upper = quantile(estimate, probs = 0.975, na.rm=T),
            estimate = mean(estimate, na.rm=T),
            sign = ifelse(upper < 0 | lower > 0, "*", "n.s.")) %>% 
  mutate(dv = factor(dv,
                     levels = c(0:4),
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))) %>% 
  
  ggplot(., aes(x = ideology, y = estimate, ymin = lower, ymax = upper, colour = treatment, lty = sign)) +
  geom_hline(yintercept = 0, lty = "dotted") +
  geom_point(position = position_dodge(width = 1), ) +
  geom_errorbar(position = position_dodge(width = 1), width = 0.1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_grid(out~dv, scales = "free_y") +
  labs(x = "Ideology", y = "Marginal effect of public (vs legislators)") +
  scale_linetype_discrete("Significance") + 
  scale_colour_discrete("Treatment", labels = c("Firms", "Worker")) + 
  scale_shape_discrete("Treatment", labels = c("Firms", "Worker")) + 
  NULL

ggsave("figures/FigureG3.pdf",
       width = 21, height = 16, units = "cm")

# . Bootstrapping for Figure G4 and G5 ####


if (!file.exists("output/bootstrapping_predictions_m2_treatment.rds")) {
  
  out_df_treatmentideology_pred <- data.frame(out = as.character(),
                                              dv = as.character(),
                                              ideology = as.numeric(),
                                              estimate = as.numeric(),
                                              id = as.numeric())
  
  nboot <- 1000
  run_start <- Sys.time()
  set.seed(123)
  
  pb = txtProgressBar(min = 0, max = nboot, initial = 0, style = 3) 
  
  for (i in 1:nboot){
    
    temp_df <- df3[sample(nrow(df3), nrow(df3), replace = T), ]
    
    m2st <- clmm2(location = subsidies ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2tt <- clmm2(location = tariffs ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2xt <- clmm2(location = tax ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    m2it <- clmm2(location = infrastructure ~ type * treatment * ideology + gender + ideology + university + treatment, random = country, Hess = T, data = temp_df)
    
    out_temp <- rbind(
      marginaleffects::predictions(m2st, by = c("subsidies", "type", "ideology", "treatment")) %>% 
        mutate(out = "Subsidies") %>% 
        rename(dv = subsidies) %>% 
        select(out, dv, type, treatment, ideology, estimate),
      marginaleffects::predictions(m2tt, by = c("tariffs", "type", "ideology", "treatment")) %>% 
        mutate(out = "Tariffs") %>% 
        rename(dv = tariffs) %>% 
        select(out, dv, type, treatment, ideology, estimate),
      marginaleffects::predictions(m2xt, by = c("tax", "type", "ideology", "treatment")) %>%
        mutate(out = "Tax cuts") %>% 
        rename(dv = tax) %>% 
        select(out, dv, type, treatment, ideology, estimate),
      marginaleffects::predictions(m2it, by = c("infrastructure", "type", "ideology", "treatment")) %>% 
        mutate(out = "Infrastructure") %>% 
        rename(dv = infrastructure) %>% 
        select(out, dv, type, treatment, ideology, estimate)) %>% 
      mutate(id = i,
             dv = as.character(dv))
    
    out_df_treatmentideology_pred <- bind_rows(out_df_treatmentideology_pred,
                                               out_temp)
    
    setTxtProgressBar(pb,i)
    
    
  }
  
  run_end <- Sys.time()
  run_end - run_start
  
  write_rds(out_df_treatmentideology_pred, "output/bootstrapping_predictions_m2_treatment.rds")
} else{
  out_df_treatmentideology_pred <- readRDS("output/bootstrapping_predictions_m2_treatment.rds")
}

# .. Figure G4 ####

out_df_treatmentideology_pred %>% 
  group_by(out, dv, type, ideology, treatment) %>% 
  summarise(lower = quantile(estimate, probs = 0.025, na.rm=T),
            upper = quantile(estimate, probs = 0.975, na.rm=T),
            estimate = mean(estimate, na.rm=T),
            sign = ifelse(upper < 0 | lower > 0, "*", "n.s.")) %>% 
  mutate(dv = factor(dv,
                     levels = c(0:4),
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))) %>% 
  filter(treatment == "workers") %>% 
  ggplot(., aes(x = ideology, y = estimate, ymin = lower, ymax = upper,shape = type)) +
  geom_hline(yintercept = 0, lty = "dotted") +
  geom_point(position = position_dodge(width = 1), ) +
  geom_errorbar(position = position_dodge(width = 1), width = 0.1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_grid(out~dv, scales = "free_y") +
  labs(x = "Ideology", y = "Predicted probability") +
  scale_colour_discrete("Treatment", labels = c("Firms", "Worker")) + 
  scale_shape_discrete("Group", labels = c("Legislator", "Public")) + 
  NULL

ggsave("figures/FigureG4.pdf",
       width = 21, height = 16, units = "cm")

# .. Figure G5 ####


out_df_treatmentideology_pred %>% 
  group_by(out, dv, type, ideology, treatment) %>% 
  summarise(lower = quantile(estimate, probs = 0.025, na.rm=T),
            upper = quantile(estimate, probs = 0.975, na.rm=T),
            estimate = mean(estimate, na.rm=T),
            sign = ifelse(upper < 0 | lower > 0, "*", "n.s.")) %>% 
  mutate(dv = factor(dv,
                     levels = c(0:4),
                     labels = c("Completely oppose", "Oppose", "Neither/nor", "Support", "Completely support"))) %>% 
  filter(treatment == "firms") %>% 
  ggplot(., aes(x = ideology, y = estimate, ymin = lower, ymax = upper,shape = type)) +
  geom_hline(yintercept = 0, lty = "dotted") +
  geom_point(position = position_dodge(width = 1), ) +
  geom_errorbar(position = position_dodge(width = 1), width = 0.1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_grid(out~dv, scales = "free_y") +
  labs(x = "Ideology", y = "Predicted probability") +
  scale_colour_discrete("Treatment", labels = c("Firms", "Worker")) + 
  scale_shape_discrete("Group", labels = c("Legislator", "Public")) + 
  NULL

ggsave("figures/FigureG5.pdf",
       width = 21, height = 16, units = "cm")

# Main analyses without control variables ####

# .. Table D1 ####


m1snc <- clmm2(location = subsidies ~ type, random = country, Hess = T, data= df3)
m1tnc <- clmm2(location = tariffs ~ type, random = country, Hess = T, data= df3)
m1xnc <- clmm2(location = tax ~ type, random = country, Hess = T, data= df3)
m1inc <- clmm2(location = infrastructure ~ type, random = country, Hess = T, data= df3)

modelsummary::modelsummary(models = list("Infrastructure" = m1inc, "Subsidies" = m1snc, "Tariffs" = m1tnc, "Tax cuts" = m1xnc),
                           output = "tables/TableD1.tex",
                           title = "Legislators, the public and policy support without controls",
                           label = "tab:model1nocontrols",
                           estimate = "{estimate}{stars}", 
                           stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
                           notes = "Note: *** p < 0.01, ** p < 0.05, * p < 0.1. Entries are unstandardised regression coefficients from a multi-level cumulative link model. Countries function as upper-level groups.",
                           fmt = 2,
                           coef_map = c('typepublic' = 'Public (ref. legislators)',
                                        '0|1' = '0|1',
                                        '1|2' = '1|2',
                                        '2|3' = '2|3',
                                        '3|4' = '3|4'),
                           gof_map = c("nobs", "aic", "bic"))

# Descriptives, samples and response rates ####

#As this data is highly sensitive, this data is not shareable. Please consult Analysis.html for the outputs.
if (dir.exists("not_shareable/")) {

population <- readRDS("not_shareable/population_parliamentarians.rds") %>% 
  filter(country %in% europe)

sample <-  readRDS("not_shareable/dat.rds") %>%
  select(contains(c("q3")), -contains("_raw"), first_name, last_name)

response_data <- dplyr::left_join(population, sample, by = c("first_name", "last_name")) %>% 
  mutate(completed = ifelse(!is.na(q3_index), 1, 0))

response_data %>% 
  select(completed, country, gender, gov_opp, parliament_chamber, age)

# .. Figure A1 ####


response_data %>% 
  mutate(age_group = cut(age,
                         breaks = seq(0, max(age, na.rm = TRUE) + 10, by = 10),
                         right = FALSE,
                         labels = paste0(seq(0, max(age, na.rm = TRUE) - 10, by = 10),
                                         "-", 
                                         seq(9, max(age, na.rm = TRUE) - 1 + 10, by = 10)))) %>% 
  group_by(age_group) %>% 
  summarise(rr = mean(completed, na.rm=T)) %>% 
  ggplot(., aes(x = age_group, y= rr)) +
  geom_bar(stat="identity", fill = "lightgray", colour = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 8)) +
  ylab("Response rate") +
  xlab("Age group") +
  NULL

ggsave("figures/FigureA1.pdf",
       width = 16, height = 9, units = "cm")

# .. Figure A2 ####


response_data %>% 
  group_by(country) %>% 
  summarise(rr = mean(completed, na.rm=T)) %>% 
  ggplot(., aes(x = country, y= rr)) +
  geom_bar(stat="identity", fill = "lightgray", colour = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 8)) +
  ylab("Response rate") +
  xlab("Country") +
  NULL

ggsave("figures/FigureA2.pdf",
       width = 16, height = 9, units = "cm")

# .. Figure A3 ####

response_data %>% 
  group_by(gender) %>% 
  summarise(rr = mean(completed, na.rm=T)) %>% 
  ggplot(., aes(x = gender, y= rr)) +
  geom_bar(stat="identity", fill = "lightgray", colour = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 8)) +
  ylab("Response rate") +
  xlab("Gender") +
  NULL

ggsave("figures/FigureA3.pdf",
       width = 16, height = 9, units = "cm")

# .. Figure A4 ####


response_data %>% 
  filter(!is.na(gov_opp)) %>% 
  group_by(gov_opp) %>% 
  summarise(rr = mean(completed, na.rm=T)) %>% 
  ggplot(., aes(x = gov_opp, y= rr)) +
  geom_bar(stat="identity", fill = "lightgray", colour = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 8)) +
  ylab("Response rate") +
  xlab("Government status") +
  NULL

ggsave("figures/FigureA4.pdf",
       width = 16, height = 9, units = "cm")

# .. Figure A5 ####

response_data %>% 
  group_by(parliament_chamber) %>% 
  summarise(rr = mean(completed, na.rm=T)) %>% 
  ggplot(., aes(x = parliament_chamber, y= rr)) +
  geom_bar(stat="identity", fill = "lightgray", colour = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 8)) +
  ylab("Response rate") +
  xlab("Chamber") +
  NULL

ggsave("figures/FigureA5.pdf",
       width = 16, height = 9, units = "cm")
} else{
  
  message("To protect respondents privacy, the data of legislators cannot be shared.")
}

